module hertzian_class
    use ellipticintegrals_module
    implicit none

    integer, parameter :: cylindricalcontact = 1
    integer, parameter :: sphericalcontact = 2
    integer, parameter :: ellipticalcontact = 3
    double precision, parameter :: pi = 3.141592653589793238462643383279502884197d0
    double precision, parameter :: errtol = 1.d-3

    type hertzian_t
        integer :: contactshape
        double precision :: Fn
        double precision :: young1, young2, poiss1, poiss2
        double precision :: rx1, rx2, ry1, ry2, alpha, width
        double precision :: effectivemodulus, rx, ry
        double precision :: ax, ay, beta, delta, p0
    end type hertzian_t




contains




function newhertzianobj()
    implicit none
    type (hertzian_t) :: newhertzianobj
    double precision :: myNaN, mm(2)
    mm = sqrt([ -1.d0, 0.d0 ])
    myNaN = mm(1)
    newhertzianobj%contactshape = 0
    newhertzianobj%Fn = myNaN
    newhertzianobj%young1 = myNaN
    newhertzianobj%young2 = myNaN
    newhertzianobj%poiss1 = myNaN
    newhertzianobj%poiss2 = myNaN
    newhertzianobj%rx1 = myNaN
    newhertzianobj%rx2 = myNaN
    newhertzianobj%ry1 = myNaN
    newhertzianobj%ry2 = myNaN
    newhertzianobj%alpha = myNaN
    newhertzianobj%width = myNaN
    newhertzianobj%effectivemodulus = myNaN
    newhertzianobj%rx = myNaN
    newhertzianobj%ry = myNaN
    newhertzianobj%ax = myNaN
    newhertzianobj%ay = myNaN
    newhertzianobj%beta = myNaN
    newhertzianobj%delta = myNaN
    newhertzianobj%p0 = myNaN
end function newhertzianobj





subroutine calchertz(c)
    implicit none
    type (hertzian_t), intent (inout) :: c
    double precision :: k, a, b, d, e, m(2,2)
    double precision :: aa, bb, valp1, valp2, l1, m1, l2, m2
    integer :: ier

    if (c%contactshape == cylindricalcontact) then
        c%effectivemodulus = 1/ ( (1-c%poiss1**2)/c%young1 + (1-c%poiss2**2)/c%young2 )
        c%rx = 2.d0 / (1/c%rx1 + 1/c%rx2)
        c%ax = sqrt(2.d0 / pi / c%width * c%rx / c%effectivemodulus)
        c%p0 = 4.d0/pi * C%Fn/c%ax/c%width
    elseif (c%contactshape == sphericalcontact) then
        c%effectivemodulus = 1/ ( (1-c%poiss1**2)/c%young1 + (1-c%poiss2**2)/c%young2 )
        c%rx = 2.d0 / (1.d0/c%rx1 + 1.d0/c%rx2)
        c%ax = (c%Fn/c%effectivemodulus*c%rx*0.375d0)**(1.d0/3.d0)
        c%p0 = 3.d0/2.d0*c%Fn/pi/c%ax/c%ax
        c%delta = 0.75d0 * c%Fn / c%effectivemodulus / c%ax
    elseif (c%contactshape == ellipticalcontact) then
        c%effectivemodulus = 1/ ( (1-c%poiss1**2)/c%young1 + (1-c%poiss2**2)/c%young2 )
        aa = 0.5d0/c%rx1 + 0.5d0/c%rx2
        bb = 0.5d0/c%ry1 + 0.5d0/c%ry2
        d = 0.5d0/c%ry2 - 0.5d0/c%rx2
        m = reshape([aa+d*sin(c%alpha)**2, -d*sin(2*c%alpha), -d*sin(2*c%alpha), bb-d*sin(c%alpha)**2 ], [2,2])
        e = sqrt( ( (m(1,1)-m(2,2)) / 2 )**2 + m(1,2)*m(2,1) )
        valp1 = (m(1,1)+m(2,2)) / 2 + e
        valp2 = valp1 - 2*e
        l1 = m(1,2)
        m1 = valp1-m(1,1)
        l2 = valp1-m(2,2)
        m2 = m(2,1)
        if ( l1*l1+m1*m1 > l2*l2+m2*m2 ) then
            c%beta = atan2(m1, l1)
        else
            c%beta = atan2(m2,l2)
        endif
        do while (c%beta > pi/2) 
            c%beta = c%beta - pi
        end do
        do while (c%beta < -pi/2) 
            c%beta = c%beta + pi
        end do
        if (c%beta < -pi/4) then
          c%beta = c%beta + pi/2
          b = valp1
          a = valp2
        elseif (c%beta > pi/4) then
          c%beta = c%beta - pi/2
          b = valp1
          a = valp2
        else
          a = valp1
          b = valp2
        end if
        c%rx = 1/a
        c%ry = 1/b
        k = contactellipticity(a/b)
        c%ax = (c%Fn /2.d0/pi / c%effectivemodulus * c%rx * RD(0.d0, 1.d0 / k / k, 1.d0, errtol, ier))**(1.d0/3.d0)
        c%ay = c%ax / k
        c%p0 = 1.5d0/pi  * c%Fn / c%ax / c%ay
        c%delta = 1.5d0/pi * c%Fn / c%effectivemodulus / c%ay * RF(k*k, 1.d0, 0.d0, errtol, ier)
    else
        print *, "Error: unknown contact shape"
        stop
    end if
end subroutine calchertz




function contactellipticity(ab)
! given the geometric pramater A/B = (1/Rx) / (1/Ry) obtain the ellipticity
! of the contact surface a/b of an elliptic contact
    implicit none
    double precision :: contactellipticity
    double precision, intent(in) :: ab
    double precision, parameter :: ktabprec = 0.001d0
    double precision, parameter, dimension(100) ::  &
    ktab = [ 1.000d0, 0.631d0, 0.483d0, 0.400d0, 0.346d0, 0.308d0, 0.279d0, 0.256d0, 0.238d0, 0.222d0, &
             0.209d0, 0.198d0, 0.189d0, 0.180d0, 0.172d0, 0.166d0, 0.160d0, 0.154d0, 0.149d0, 0.144d0, &
             0.140d0, 0.136d0, 0.132d0, 0.129d0, 0.126d0, 0.123d0, 0.120d0, 0.118d0, 0.115d0, 0.113d0, &
             0.110d0, 0.108d0, 0.106d0, 0.104d0, 0.103d0, 0.101d0, 0.099d0, 0.098d0, 0.096d0, 0.095d0, &
             0.093d0, 0.092d0, 0.091d0, 0.089d0, 0.088d0, 0.087d0, 0.086d0, 0.085d0, 0.084d0, 0.083d0, &
             0.082d0, 0.081d0, 0.080d0, 0.079d0, 0.078d0, 0.077d0, 0.077d0, 0.076d0, 0.075d0, 0.074d0, &
             0.074d0, 0.073d0, 0.072d0, 0.072d0, 0.071d0, 0.070d0, 0.070d0, 0.069d0, 0.068d0, 0.068d0, &
             0.067d0, 0.067d0, 0.066d0, 0.066d0, 0.065d0, 0.065d0, 0.064d0, 0.064d0, 0.063d0, 0.063d0, &
             0.062d0, 0.062d0, 0.061d0, 0.061d0, 0.060d0, 0.060d0, 0.060d0, 0.059d0, 0.059d0, 0.058d0, &
             0.058d0, 0.058d0, 0.057d0, 0.057d0, 0.057d0, 0.056d0, 0.056d0, 0.056d0, 0.055d0, 0.055d0 ]
    double precision :: rab, high, low, mid, midresid
    integer :: ier
    ! normalize
    rab = abs(ab)
    rab = max(rab, 1.d0/rab)
    if (rab < 1) then
        print *, "rab < 1"
        stop
    end if
    ! bisection method
    if (ceiling(rab) <= size(ktab)) then
        high = ktab(floor(rab)) + ktabprec
        low = ktab(ceiling(rab)) - ktabprec
    else
        high = ktab(size(ktab)) + ktabprec
        low = 0.d0
    end if
    mid = low + (high-low)/2.d0
    do while (mid /= low .and. mid /= high)
        midresid = rab - RD(0.d0, 1.d0, mid*mid, errtol, ier) / RD(0.d0, mid*mid, 1.d0, errtol, ier)
        if (midresid > 0.d0) then
            high = mid
        else
            low = mid
        end if
        mid = low + (high-low)/2.d0
    end do
    ! de-normalize
    if (rab > abs(ab)) then
        mid = 1.d0/mid
    end if
    if (ab < 0.d0) then
        mid = -mid
    end if
    contactellipticity = mid
    return
end function contactellipticity




end module hertzian_class
